clear all; clc; close all;
RC = 1e-6;
wmax=1/RC*2*pi*2;
dw=wmax/100;
w=0:dw:wmax;

a = [RC 1];  %[a1 a0]
b = [1];     %[b0]

figure(1)   
freqs(b, a, w);  %вывести АЧХ/ФЧХ, приняв по оси частот w
H = freqs(b,a,w);% посчитать ПФ, но график не выводить

figure(2)
subplot(2,1,1)
plot(w/1/pi/1e6, abs(H)); xlabel('f, MHz'); ylabel('|H|');
subplot(2,1,2);
plot(w/1/pi/1e6, unwrap(angle(H))); xlabel('f, MHz'); ylabel('arg H');

%Находим импульсную характеристику
sys=tf(b,a);
[y,t]=impulse(sys);
figure(3)
plot(t,y); xlabel('t,sec'); ylabel('h(t)');

%Находим цифровой аналог
Fs = 2*wmax/2/pi;
[bz, az] = bilinear(b, a, Fs);
[Hz, wz] = freqz(bz, az);
figure;
plot(w/2/pi/1e6, abs(H), Fs*wz/2/pi/1e6, abs(Hz));
xlabel('f, MHz'); ylabel('|H|')

% Шум
Td = 1/Fs;
Tmod = 1000*Td;
t = 0:Td:Tmod;
stdn = 13; n = stdn*randn(1, length(t));
y = filter(bz, az, n);
figure;
plot(t, [n; y]) % plot(t, x, t, y)
xlabel('t, sec');


nf = fft(n);
yf = fft(y);
mn = sqrt(mean(abs(nf).^2))*2.5;
nf = nf / mn;
yf = yf / mn;
f = 0: 1/Tmod:(1/Td);
figure;
plot(f/1e6, [abs(nf); abs(yf)]/max(abs(yf)), w/2/pi/1e6, abs(H), 'LineWidth', 3);
xlim([0 Fs/1e6/2])
xlabel('f, Hz')